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Abstract 


Genotype by environment interaction and stability analyses are among the 

most important evaluations conducted in plant breeding. In this study, we 
evaluated the yield-related traits of 32 Tunisian barley (Hordeum vulgare L.) 
accessions over three consecutive cropping seasons in a semi-arid environment. 
Phenotypic analysis identified heading date and spike length as the two major 
traits contributing most to the total phenotypic variation under a semiarid 
climate. Hierarchical clustering grouped the 32 accessions into four groups. 
Although the effect of the interaction between genotype and environment was 
important for yield (48%), it had comparatively little influence on heading date 
(9.9%) and plant height (8.14%). Stability analysis identified the djebali accession, 
with the smallest coefficient of variability, as a stable genotype across the three 
assessed environments. Overall, based on the higher yield and small coefficient of 
variance, we selected 11 promising genotypes. In addition, varieties developed by 
the Tunisian breeding program were found to show high mean yield, stability 
across all environments, and greater adaptability. Accessions with superior 
adaptation and stability will be introduced into the national breeding program 
for further evaluation and characterization. 
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1. Introduction 


Barley (Hordeum vulgare L.) is one of the oldest cultivated crops worldwide, ranking 
fourth after wheat, rice, and maize. In Tunisia, barley is mainly cultivated in regions 
with arid and semiarid climates that receive less than 400 mm of rainfall annually, 
and as part of the Tunisian breeding program for barley, genetic studies are currently 
being conducted with the aim of enhancing yields. Yield in barley is a complex trait 
governed by several genes that interact with the environment, and consequently, 

the selection of genotypes based on performance in a single environment is an 
ineffective approach for varietal selection (Shrestha et al., 2012). 


In this regard, the analysis of genotype by environment interactions is of particular 
importance, notably in regions of North Africa, wherein barley is often cultivated 
under adverse conditions of drought, high temperatures, and irregular rainfall 

(van Oosterom & Ceccarelli, 1993). Accordingly, trials encompassing multiple 
environments are required to identify optimal environments for selecting genotypes 
for enhanced grain yield (Gauch & Zobel, 1997). In this respect, the most widely 
adopted approaches in the study of genotype by environment interactions are based 
on models incorporating genotype main effects and Genotype x Environment 
interaction effects (Yan et al., 2000). 
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In this study, we analyzed genotype by environment interactions for a collection 

of 32 Tunisian barley accessions that were grown over three consecutive growing 
seasons in a single location. In order to prioritize the component traits important for 
further genetic dissection and improvement, we assessed the phenotypic diversity of 
important yield-related traits. 


2. Material and Methods 
2.1. Germplasm and Phenotyping 


For the purposes of this study, we assessed 32 Tunisian barley lines (Hordeum 
vulgare L.) comprising four cultivars developed by the Tunisian breeding program 
in collaboration with ICARDA (Rihane, Manel, Lemsi, and Kounouz), one accession 
with uncertain improvement status, and 27 landraces. All accessions were obtained 
from the U.S. National Plant Germplasm System international database. According 
to passport data, the 28 accessions other than the aforementioned four cultivars were 
collected or donated from Tunisia between 1922 and 1972 (Table 1). During the 
three cropping seasons of 2016/2017, 2017/2018, and 2018/2019, field trials were 
conducted in El Kef, Tunisia, which is characterized by a semiarid climate with an 
average annual rainfall of 380 mm (Table 2). Each accession was planted in two 
replicate plots, each comprising two 2.5-m-long rows planted with 65 seeds with 

an inter-row spacing of 25 cm. With the exception of yield, all assessed traits were 
measured from five plants. The traits evaluated in this study were as follows: grain 
yield: the weight of grain harvested from the two rows; plant height: measured from 
the soil surface to tip of the spike (excluding awns); days to heading: the number 

of days from sowing to the time when 50% of the ear had emerged; peduncle 

length: measured from the base of the spike to the flag leaf node; internode length: 
measured from the flag leaf node to the next node; awn length: the length of the awn 
in the central spikelet; spike length: measured from the base of spike to the tip of 
the terminal spikelet (excluding awns); and kernel number per spike: the number of 
grains per spike. 


2.2. Statistical Analysis of Phenotypic Data 


For each trait, descriptive statistical data were obtained based on the average barley 
accession data. We used the following additive main effects and multiplicative 
interaction model that combines two standard methods, namely, analysis of variance 
(ANOVA) with principal component analysis (PCA) analysis (Zobel et al., 1988): 


n 
Yijr =W+ Gi + Gj + by(ej) + Dy AKHiKYjk + Pi + Ey: 
k=1 
where Yj, is the phenotypic trait, 1 is the grand mean, g; is the main effect of the 
ith genotype (G), and e; is the main effect of the jth environment (E). A, is the 
singular value for the interaction principal component (IPC) axis, k, ox and y;x are 
the genotype and environment IPC scores (i.e., the left and right singular vectors) 
for axis; k, b,(e;) is the effect of block r within environment j, r is the number of 
blocks (in this study, the number of blocks was not considered, and thus, the effect 
of blocks was removed), pj is the residual containing all multiplicative terms not 
included in the model, n is the number of axes or principal components retained 
by the model, and ¢;, is the experimental error, assumed to be independent of the 
distribution ej ~ N (0, o). 
Analysis of the additive main effects and multiplicative interaction model and the 
genotype main effects and Genotype X Environment interaction effects model were 
performed using Genotype x Environment Analysis with R software in Windows 
(Pacheco et al., 2016). 


3. Results 
3.1. Descriptive Statistics 


Basic statistical parameters for each of the three growing seasons (mean, standard 
deviation, minimum, and maximum values) for the barley accessions are 
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Table 1 A list of the 32 Tunisian barley accessions used for phenotyping. 








Number Accession name Status Origin 
Gl Rihane Cultivar Tunisia 
G2 Kounouz Cultivar Tunisia 
G3 Lemsi Cultivar Tunisia 
G4 Manel Cultivar Tunisia 
G5 175 Uncertain Ariana 
G6 2528-23 Landrace Siliana 
G7 3124-8 Landrace Siliana 
G8 djebali Landrace Manouba 
G9 djebali Landrace Manouba 
G10 djebali Landrace Manouba 
Gll djebali Landrace Manouba 
G12 frigui Landrace Kebili 
G13 frigui Landrace Kebili 
G14 djebali Landrace Kebili 
G15 1110-30 Landrace Kebili 
G16 jebali Landrace Kebili 
G17 jebali Landrace Kebili 
G18 djebali Landrace Kebili 
G19 hmira Landrace Kebili 
G20 djebali Landrace Kebili 
G21 jebali Landrace Kebili 
G22 djebali Landrace Kebili 
G23 frigui Landrace Kebili 
G24 jebali Landrace Kebili 
G25 jebali Landrace Kebili 
G26 djebali Landrace Kebili 
G27 jebali Landrace Bizerte 
G28 tounsi Landrace Tozeur 
G29 safra Landrace Tozeur 
G30 commune A Landrace Unknown 
G31 cowra Landrace Unknown 
G32 cowra Landrace Unknown 





Table 2 Descriptions of the environmental conditions during three evaluated cropping seasons. 











Years Nov. Dec. Jan. Feb. Mar. Apr. May 
Average T (°C) Vegetative stage Flowering stage _—_ Grain filling 
2018/2019 12.6 12.5 9.0 7.0 7.3 10.6 13.4 15.9 
2017/2018 12.2 11.0 7.9 9.7 8.1 12.2 15.9 25.0 
2016/2017 12.6 13.4 10.5 6.3 9.4 11.9 15.7 20.8 
Precipitation(mm) Vegetative stage Flowering stage _—_ Grain filling 
2018/2019 530.0 48.3 66.8 115.0 74.0 80.0 47.0 99.3 
2017/2018 309.0 82.5 36.8 32.5 31.7 41.2 36.3 48.5 
2016/2017 181.3 44.2 53.1 35.0 24.9 4.3 17.5 2.3 
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Table 3 Descriptive statistics of the studied traits. 








Traits Years Minimum Maximum Mean SD 
HD (heading date) 2017 118.0 135.0 123.7 4.7 
2018 122.0 133.0 126.4 3.6 
2019 125.0 141.0 131.8 3.9 
PH (plant height) 2017 55.0 83.0 68.2 6.5 
2018 72.0 90.0 80.7 3.9 
2019 92.0 114.0 102.5 5.4 
INL (internode length) 2017 4.0 16.0 8.8 2.8 
2018 6.0 19.0 10.6 3.5 
2019 12.0 24.0 18.3 2.8 
PL (peduncle length) 2017 6.0 18.0 12.3 2.5 
2018 8.0 19.0 14.3 2.6 
2019 17.6 34.0 26.3 3.9 
SL (spike length) 2017 5.0 9.0 7.6 11 
2018 5.0 9.0 7.5 1.0 
2019 5.0 10.0 7.9 0.9 
KNS (kernel number per spike) 2017 52.0 72.0 62.4 5.6 
2018 54.0 78.0 65.3 5.7 
2019 64.0 78.0 69.6 4.4 
AwL (awn length) 2017 1.0 14.0 10.7 2.2 
2018 1.0 15.0 10.1 22 
2019 1.0 16.0 12.2 2.6 
YLD (yield) 2017 174.0 526.0 368.5 81.4 
2018 128.0 660.0 394.6 109.8 
2019 340.0 623.0 457.8 62.5 





Table 4 Pearson coefficients of correlation between studied traits. 








Traits HD PH INL PL SL KNS AwL 
PH 0.10 

INL 0.04 0.15 

PL —0.54 0.18 0.31 

SL —0.52* —0.23 —0.26 0.29 

KNS —0.16* 0.10 0.28 0.03 0.16 

AwL —0.39* 0.22 —0.26 0.34 0.41* -0.15 

YLD —0.26 0.45* 0.06 0.18 —0.08 0.27 0.17 





* The correlation is significant at the p < 0.05. HD — heading date; PH — plant height; INL — internode length; PL - peduncle length; SL — spike length; KNS 
- kernel number per spike; AwL - awn length; YLD - yield. 


summarized in Table 4. Almost all traits showed an increase in mean values over 
the three cropping seasons; for example, mean grain yield increased from 368.5 g 
in 2017 to 457.8 g in 2019, and plant height increased from 68 cm to 102.5 cm. In 
contrast, there was a delay in the heading date from 123 days in 2017 to 131 days in 
2019. 


PCA analysis revealed that 68.7% of the phenotypic variation could be explained 

by the first three components, with the first two components, PCA1 and PCA2, 
explaining 28.7% and 24.9% of the total variation, respectively (Figure 1A). Heading 
date and spike length were found to be major contributors to these two components, 
whereas among the assessed accessions, kernel number per spike exhibited the 
lowest contribution to the total phenotypic variation. Hierarchical clustering based 
on the mean value of the target traits over the three cropping seasons grouped the 
accessions into four clusters (Figure 1B). The Lemsi cultivar was considered to 
represent a discrete group, on account of its late heading and awnless spike. The 
other modern varieties were grouped in the second cluster, the constituents of 
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Figure 1 Principal component (A) and hierarchical cluster (B) analyses. The first two 
principal components, PCA1 and PCA2, explained 28.7% and 24.9% of the total variation. 
Days to heading (DH) and spike length (SL) make the major contributions to the total 
variation. Hierarchical cluster based on the mean values of the studied traits across the 
three cropping seasons identified clusters of accessions. See the text for explanation of 
other abbreviations. 


which were characterized by tall plants with higher yields. The third cluster group 
comprised medium-tall plants with low yield and late heading, whereas the fourth 
group contained early heading short-statured accessions. 


The phenotypic correlations of the studied traits across the three cropping seasons 
are presented in Table 4. The combined data obtained over the three seasons 
revealed that heading date was negatively correlated with almost all other traits, with 
the exception of plant height. Plant height was found to be significantly positively 
correlated with grain yield (r = 0.45, p < 0.05), but it was negatively correlated with 
spike length. 


3.2. Analysis of Variance 


The results of analyses of variance for the grain yield-related traits of each genotype 
in each growth season based on interaction effects are presented in Table 5. Total 
variance among the genotypes, environment (seasons), and interaction effects 

was found to be highly significant (p < 0.01). For heading date, spike length, 
kernel number per spike, awn length, and yield, genotypic effects were found 

to be considerably stronger than the environmental effects, ranging from 35.6% 
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Table 5 Analysis of variance across three cropping seasons. 








Traits Genotype variance Environment variance Genotype/Environment Significance (p < 
(G) (E) variance (G x E) 0.01) 
Heading date (HD) 48.0 42.0 9.9 0 
Plant height (PH) 3.7 88.0 8.1 0 
Internode length (INL) 20.8 64.7 14.4 0 
Peduncle length (PL) 6.1 81.3 12.5 0 
Spike length (SL) 66.9 27 30.4 0 
Kernel number per spike (KNS) 42.1 24.6 33.2 0 
Awn length (AwL) 65.0 14.0 20.4 0 
Grain Yield (Yld) 35.6 16.3 48.0 0 





to 66.9%. The effect of the interaction between genotype and environment on 

yield was important for yield (48%), although similar effects were relatively small 
with respect to heading date (9.9%) and plant height (8.14%). A large phenotypic 
variation explained by genotypes indicated the diversity of the genotypes and that a 
major fraction of the observed variation in heading date and spike length could be 
attributed to genetic effects. 


3.3. Analysis of Stability and Genotype by Environment Interaction Effects 


Stability parameters are useful for characterizing genotypes by indicating their 
relative performance in different environments. In the present study, we combined 
the data for yield and yield-related traits for all three cropping seasons with stability 
statistics by measuring the Shukla variance (Shukla et al., 1972) (Table 6). 


It was accordingly found that accession G20 with the smallest coefficient of 
variability (-207) was the most stable genotype in the selection ranks for yield across 
the three cropping seasons; however, there were a further 11 genotypes selected 
based on higher yield and small coefficient of variance (Rihane, Manel, G12, G13, 
G14, G16, G23, G25, G29, G31, and G32). 


The number of genotypes selected based on coefficient of variability (s*4;) ranking 
varied among the traits. Two accessions (G13 and G17) were found to be stable with 
respect to heading date; two (G18 and G20) were stable for plant height; four were 
stable for spike length (G7, G24, G25, and G30); and two were stable in terms of awn 
length (G3 and G27). 


The classification of genotype instability attributable to environmental variability 
highlighted the need for further detailed studies on the behavior of these accessions 
during selection of the best genotypes. To visualize the performance of different 
genotypes in a given environment, we constructed biplots. The first two main factors 
(components) of the biplot analysis was applied to the Genotype x Environment 
interaction and explained 83.9% of the total variation, thereby indicating that 
biplot graphics explain most of the sums of squares and Genotype X Environment 
interactions with respect to genotype. The lines (blue lines) that connect the test 
environments to the biplot origin are referred to as environment vectors. The cosine 
of the angle between the vectors of the two environments approximates to the 
correlation between these environments. Environment 2017 (E2017) and E2018 
were found to be positively correlated (acute angle) and E2018 and E2019 were 
highly correlated, whereas E2017 and E2019 appeared to be uncorrelated (right 
angle) (Figure 2). The close associations detected among test environments indicate 
that we can obtain similar information regarding genotypes. The concentric circles 
on the biplot enable visualization of the length of the environmental vectors, which 
provides a measure of environmental discrimination. Accordingly, the two years 
2017 and 2018 show similar vector lengths and are thus considered discriminating 
and informative environments, whereas 2019 was found to be less informative. The 
smaller the angle of the environmental vector with the AED (average-environment 
axis) (Figure 2), the more representative it is of other test environments. Thus, 
E2017 and E2018 were found to be the most representative, whereas E2019 was 
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Table 6 Analysis of stability based on measurements of regression (bj) and variability (s? 4;) coefficients. 








GEN Mean (g) SD CV (%) sai R2 R? 

G20 370.00 11.65 3.14 —207.34 0.99 3,322.29 
G18 398.33 16.92 4.24 -196.38 0.98 3,984.45 
G12 465.16 36.72 7.89 —192.98 0.99 7,064.26 
G30 285.16 112.71 39.52 -173.06 0.99 4,552.37 
G16 480.83 33.00 6.86 -80.31 0.94 6,335.50 
G7 304.00 106.98 35.19 108.24 0.98 3,836.80 
G22 402.16 55.34 13.76 134.50 0.94 35.04 

G5 443.83 94.60 2131 180.25 0.97 2,415.72 
G29 462.16 80.75 17.47 555.46 0.94 1,315.56 
G17 378.16 25.79 6.82 911.43 0.15 3,755.50 
G19 401.66 41.59 10.35 1,122.32 0.61 7,079.35 
Gl 501.16 89.51 17.86 1,289.14 0.90 2,232.80 
G3 378.58 101.47 26.80 1,401.54 0.92 3,474.08 
G14 466.50 33.72 7.22 2,034.06 0.01 2,852.25 
G4 414.63 63.60 15.33 2,247.34 0.69 1,150.69 
G2 416.33 88.63 21.28 3,011.98 0.79 2,673.23 
G25 418.83 41.06 9.80 3,107.45 0.01 4,363.79 
G9 334,33 114.50 34.25 3,401.61 0.86 5,604.71 
G3l 452.16 79.84 17.65 3,579.41 0.70 2,277.93 
G32 434.50 53.20 12.24 4,145.73 0.23 2,550.20 
G28 350.33 119.58 34.13 4,836.63 0.82 6,657.84 
G13 439.16 57.69 13.13 6,358.17 0.01 4,921.38 
G23 429.16 77.92 18.15 6,556.40 0.44 3,431.88 
Gll 436.16 136.09 31.20 6,774.87 0.81 9,781.99 
G6 341.83 92.28 26.99 7,427.06 0.55 4,404.46 
G27 365.00 70.70 19.37 7,715.38 0.20 4,212.93 
G24 367.00 97.61 26.59 12,639.61 0.32 6,741.05 
G15 382.50 97.89 25.59 18,686.74 0.01 13,400.25 
G21 411.83 101.18 24.56 19,177.01 0.05 10,668.70 
G10 316.16 174.06 55.05 28,025.82 0.53 21,893.03 
G8 503.16 143.43 28.50 33,799.63 0.17 18,127.48 
G26 471.33 162.94 34.57 52,773.19 0.00 29,596.84 





bj - coefficient of regression; s? 4 — coefficient of variability; CV - coefficient of variance; R2 , R;2- determination coefficients. 


the least representative. Environments such as E2017 and E2018 that are both 
discriminating and informative are considered good environments for selecting 
generally adapted genotypes. The pattern of the environments in the aforementioned 
biplots is indicative of two mega-environments, one of which one is a semiarid 
climate represented by E2017 and E2018, and the other is a subhumid climate 
represented by E2019. 


A set of perpendicular lines divided the which-won-where biplot into multiple 
groups (Figure 3). A polygon was initially drawn to encompass all genotypes. 
Thereafter, lines perpendicular to each side of the polygon were drawn, starting 
from the biplot origin. The perpendicular lines radiating from the origin of the 
biplot divide the plot area as well as the test location into sectors. Note the red 

lines and genotypes in the vertex of the polygon, and the lines dividing the sectors, 
indicating the mega-environments. The best genotype in each mega-environment 
corresponds to the vertex; for example, in the sector representing E2018 and E2019, 
the genotypes G8 and G26 have the best yield, whereas in the sector representing 
E2017, G12 and G16 have the best yields. 


We subsequently evaluated the 32 assessed genotypes for mean performance and 
stability across environments (Figure 4). The single-arrowed line is the average- 
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Figure 2 Discriminativeness vs. representativeness of the three evaluated cropping 
seasons. The first two main factors (components) of the biplot analysis applied to 
Genotype X Environment interaction effects explained 83.9% of the total variation. 
Environments 2017 and 2018 were positively correlated (acute angle) and 2018 and 2019 
were highly correlated, whereas 2017 and 2019 were uncorrelated (right angle). 
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Figure 3 A which-won-where biplot. In the sector within which environments 2018 and 
2019 are plotted, the genotypes G8 and G26 have the highest yields, whereas in the sector 
in which environment 2017 is plotted, G12 and G16 have the highest yields. 


environment coordination abscissa, which indicates a higher mean yield across 
environments. Accordingly, the plot shows that accessions G8 and G26 had the 
highest mean yield, followed by G12 and G16, whereas G10 and G30 had the lowest 
mean yields. 


The higher the genotype projection on the axis, the lower is the genotype by 
environment interaction, and consequently, the more stable the genotype. Thus, 
genotypes G8 and G26 were identified as being highly unstable, whereas G20 and 
G22 are considered highly stable. 


4. Discussion 


This study was conducted at the El Kef Research Station, which is located in a region 
characterized by low and irregular rainfall throughout the barley cropping season, 
particularly at the grain filling stage, which affects grain yield, as shown in Table 2. 
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Figure 4 A Genotype x Environment interaction effects biplot representing the Means x 
Stabilities that indicate the yield rankings of the 32 assessed barley accessions. Accessions 
G8 and G26 have the highest mean yields, followed by G12 and G16, whereas G10 and 
G30 have the lowest mean yield. Genotype G8 and G26 were highly unstable, whereas G20 
and G22 were highly stable. 


During the period of flowering and grain filling, the total rainfall increased from 
20 mm in 2017 to 84.4 mm in 2018 and 146.3 mm in 2019. Given this irregular 
distribution of rainfall, we considered each year to represent a separate environment. 


Principal component analysis (PCA) was performed based on seven selected yield- 
related traits. Analysis results indicated three eigenvalues greater than 1. The first 
two components, PCA1 and PCA2, explained 28.7% and 24.9% of the total variation, 
respectively, to which heading date and spike length made the major contributions. 
Heading date was found to be negatively associated with grain yield, indicating that 
barley cultivars with a long period prior to heading are unsuitable for cultivating 

in the semiarid conditions of the El Kef region. Grain filling of late heading date 
cultivars often occurs in such unfavorable environments characterized by high 
temperatures and water deficit. In this regard, Przulj and Moméilovi¢ (2006) have 
stated that barley cultivars with very early heading have lower yield potentials. 
Accordingly, cultivars with moderately early heading tend to be characterized by 
optimal phenological development and are more suitable for large-scale production 
(Mirosavljevic¢ et al., 2016). 


Analysis of variance for grain yield-related traits revealed that the genotypic effect 
was considerably stronger than that of the environment for heading date and 
yield, with a significant Genotype x Environment interaction effect for grain yield 
(Schmalenbach et al., 2009; von Korff et al., 2008). 


Collectively, the findings of this study reveal the importance of interaction and 
stability analyses with respect to the evaluation of genotypic yield potential. By 
evaluating a range of barley accessions over three consecutive cropping seasons, 
we were able to select genotype G20 as the best accession in terms of yield stability 
across environments. In general, varieties developed by the Tunisian breeding 
program showed a high mean yield, stability across all environments, and greater 
adaptability. Selected accessions will be introduced into the Tunisian breeding 
program for further evaluation and characterization. 
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